package scu.maqiang.numeric;

import java.util.Arrays;

public class EigenSSolver extends Solver{

	/**
	 * 默认构造函数
	 */
	public EigenSSolver() {
	}

	/**
	 * 构造函数，设置求解器容许精度与最大迭代次数
	 * 
	 * @param eps          容许精度
	 * @param maxIteration 最大迭代次数
	 */
	public EigenSSolver(double eps, int maxIteration) {
		this.eps = eps;
		this.maxIteration = maxIteration;
	}

	/**
	 * 用带平移的逆迭代次求解对称特征值问题, 其中求解器使用预条件CG法 $\fAx = \lambda M x\f$
	 * @param A
	 * @param M
	 * @param x0
	 * @param alpha
	 * @param num
	 * @param D
	 * @param V
	 */
	public void inversePowerTran_PCG(SRMatrix A, SRMatrix M, double[] x0, double alpha, int num, double[] D,
			double[][] V) {
		int i, j, k;
		int n = A.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		SRMatrix L = A.add(alpha, M);
		NewIterSSolver solver = new NewIterSSolver(L);
		for (i = 0; i < num; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = x0[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					double temp = MVO.dot_product(V[j], M.mul(x));
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += V[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				double xMx = Math.sqrt(MVO.dot_product(x, M.mul(x)));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				y = M.mul(x);
				solver.PCGSSOR(y, xx, 1.5, 0);
				System.out.println("Find\t" + (i + 1) + "\t of\t" + num + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				MM = MVO.dot_product(xx, M.mul(xx));
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println(i + "  of  " + num + "  finished!");
			// write(*, *) rr - alpha
			D[i] = rr - alpha;
			for (int kk = 0; kk < n; kk++) {
				V[i][kk] = xx[kk] / Math.sqrt(MM);
			}
		}
	}

	/**
	 * 
	 * @param A
	 * @param M
	 * @param nBand
	 * @param x0
	 * @param alpha
	 * @param num
	 * @param D
	 * @param V
	 */
	public void inversePowerTran_PLU(SRMatrix A, SRMatrix M, int nBand, double[] x0, double alpha, int num, double[] D,
			double[][] V) {
		int i, j, k;
		int n = A.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		SRMatrix L = A.add(alpha, M);
		DirectSSolver solver = new DirectSSolver(L);
		for (i = 0; i < num; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = x0[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					double temp = MVO.dot_product(V[j], M.mul(x));
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += V[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				double xMx = Math.sqrt(MVO.dot_product(x, M.mul(x)));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				y = M.mul(x);
				xx = solver.LUSolve(y, 1, 0.5);
				System.out.println("Find\t" + (i + 1) + "\t of\t" + num + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				MM = MVO.dot_product(xx, M.mul(xx));
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println((i + 1) + "  of  " + num + "  finished!");
			// write(*, *) rr - alpha
			D[i] = rr - alpha;
			for (int kk = 0; kk < n; kk++) {
				V[i][kk] = xx[kk] / Math.sqrt(MM);
			}
		}
	}
	
	/**
	 * 
	 * @param A
	 * @param M
	 * @param x0
	 * @param alpha
	 * @param num
	 * @param D
	 * @param V
	 */
	public void inversePowerTran_GMRES(SRMatrix A, SRMatrix M, double[] x0, double alpha, int num, double[] D,
			double[][] V) {
		int i, j, k;
		int n = A.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		SRMatrix L = A.add(alpha, M);
		NewIterSSolver solver = new NewIterSSolver(L);
		solver.setMaxIteration(8000000);
		for (i = 0; i < num; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = x0[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					double temp = MVO.dot_product(V[j], M.mul(x));
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += V[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				double xMx = Math.sqrt(MVO.dot_product(x, M.mul(x)));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				y = M.mul(x);
				// solver.GMRES_SSOR(y, xx, 1.5, 200, 2);
				// solver.BiCGSTAB(y, xx);
				solver.CGNR(y, xx, 0);
				System.out.println("Find\t" + (i + 1) + "\t of\t" + num + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				MM = MVO.dot_product(xx, M.mul(xx));
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println((i + 1) + "  of  " + num + "  finished!");
			// write(*, *) rr - alpha
			D[i] = rr - alpha;
			for (int kk = 0; kk < n; kk++) {
				V[i][kk] = xx[kk] / Math.sqrt(MM);
			}
		}
	}

	/**
	 * 
	 * @param A
	 * @param M
	 * @param x0
	 * @param alpha
	 * @param num
	 * @param D
	 * @param V
	 */
	public void inversePowerTran_PCGSSOR(SRMatrix A, SRMatrix M, double[] x0, double alpha, int num, double[] D, double[][] V) {
		int i, j, k;
		int n = A.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		double temp = 0.0;
		SRMatrix L = A.add(alpha, M);
		NewIterSSolver solver = new NewIterSSolver(L, true);
		for (i = 0; i < num; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = x0[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					M.mul(x, y);
					temp = MVO.dot_product(V[j], y);
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += V[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				M.mul(x, y);
				double xMx = Math.sqrt(MVO.dot_product(x, y));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				M.mul(x, y);
				solver.PCGSSOR(y, xx, 1.5, 0);
				System.gc();
				System.out.println("Find\t" + (i + 1) + "\t of\t" + num + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				M.mul(xx, y);
				MM = MVO.dot_product(xx, y);
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println((i + 1) + "  of  " + num + "  finished!");
			// write(*, *) rr - alpha
			D[i] = rr - alpha;
			for (int kk = 0; kk < n; kk++) {
				V[i][kk] = xx[kk] / Math.sqrt(MM);
			}
			System.gc();
		}
	}
	
	public void subspaceIteration(SRMatrix A, SRMatrix M, double[] x0, int num, double[] D, double[][] V) {
		int n = A.getM();
		double[][] x = MVO.rand(num, n);
		//double[][] x = {{1, 1, 1, 1}, {1, 1, 1, 0.5}};
		double[][] y = new double[num][];
		for(int i = 0; i < num; i++) {
			y[i] = M.mul(x[i]);
		}
		int k = 1;
		double[][] xx = new double[num][n];
		double[][] y1 = new double[num][];
		NewIterSSolver solverA = new NewIterSSolver(A);
		solverA.setEps(1.0e-16);
		double[][] KK = new double[num][num];
		double[][] MM = new double[num][num];
		double[] xx0 = MVO.rand(num);
		EigenDSolver edsolver = new EigenDSolver();
		double[] DDi = new double[num];
		double[][] Vi = new double[num][num];
		boolean br = false;
		while(k < maxIteration) {
			for(int i = 0; i < num; i++) {
				solverA.PCGSSOR(y[i], xx[i], 1.5, 0);
			}
			MVO.matmul(xx, MVO.transpose(y), KK);
			for(int i = 0; i < num; i++) {
				y1[i] = M.mul(xx[i]);
			}
			MVO.matmul(xx, MVO.transpose(y1), MM);
			edsolver.inversePowerTran(KK, MM, xx0, 0.0, num, D, Vi);
//			Vi[0] = new double[]{-1.000000000000000, 0.506142655577088};
//			Vi[1] = new double[]{-0.751545155250783, 1.000000000000000};
			if (k > 1) {
				br = false;
				for(int i = 0; i < num; i++) {
					if(Math.abs(DDi[i] - D[i]) > 1.0e-20) {
						br = true;
						break;
					}
				}
			}
			else {
				br = true;
			}
			if(!br) {
				break;
			}
			k++;
			for(int i = 0; i < num; i++) {
				y[i] = MVO.matmul(MVO.transpose(y1), Vi[i]);
				DDi[i] = D[i];
			}
			//System.out.println(Arrays.toString(D));
		}
		//System.out.println(Arrays.toString(D));
		double[][] VVi = new double[num][];
//		System.out.println(MVO.toString(xx));
		for(int i = 0; i < num; i++) {
			VVi[i] = MVO.matmul(MVO.transpose(xx), Vi[i]);
		}
		for(int i = 0; i < num; i++) {
			double vMv = MVO.dot_product(VVi[i], M.mul(VVi[i]));
			for(int j = 0; j < n; j++) {
				V[i][j] = VVi[i][j] / Math.sqrt(vMv);
			}
		}
		
	}

	public void inversePowerTran_DSS(SRMatrix lhsMat, SRMatrix rhsMat, double[] initVal, double offsetVal, int numOfEV,
									 double[] eigenValues, double[][] eigenVectors) {
		int i, j, k;
		int n = lhsMat.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		double temp = 0.0;
		SRMatrix L = lhsMat.add(offsetVal, rhsMat);
		DirectSSolver solver = new DirectSSolver(L);
		solver.Factorize_Cholesky(1);
//		NewIterSSolver solver = new NewIterSSolver(L, true);
		for (i = 0; i < numOfEV; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = initVal[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					rhsMat.mul(x, y);
					temp = MVO.dot_product(eigenVectors[j], y);
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += eigenVectors[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				rhsMat.mul(x, y);
				double xMx = Math.sqrt(MVO.dot_product(x, y));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				rhsMat.mul(x, y);
				//solver.PCGSSOR(y, xx, 1.5, 0);
				xx = solver.Solve_Cholesky(y, 0);
				System.out.println("Find\t" + (i + 1) + "\t of\t" + numOfEV + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				rhsMat.mul(xx, y);
				MM = MVO.dot_product(xx, y);
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println((i + 1) + "  of  " + numOfEV + "  finished!");
			// write(*, *) rr - alpha
			eigenValues[i] = rr - offsetVal;
			for (int kk = 0; kk < n; kk++) {
				eigenVectors[i][kk] = xx[kk] / Math.sqrt(MM);
			}
			System.gc();
		}
	}

	public void inversePowerTran_DSS_2(SRMatrix lhsMat, SRMatrix rhsMat, double[] initVal, double offsetVal, int numOfEV,
									 double[] eigenValues, double[][] eigenVectors) {
		int i, j, k;
		int n = lhsMat.getM();
		double[] x = new double[n];
		double[] xx = new double[n];
		double[] y = new double[n];
		double[] summ = new double[n];
		double r = 0.0, rr = 0.0;
		double KK = 0.0, MM = 0.0;
		double temp = 0.0;
		SRMatrix L = lhsMat.add(offsetVal, rhsMat);
		DirectSSolver solver = new DirectSSolver(L);
		solver.Factorize_LU(1, 1.0e-10);
//		NewIterSSolver solver = new NewIterSSolver(L, true);
		for (i = 0; i < numOfEV; i++) {
			for (int kk = 0; kk < n; kk++) {
				x[kk] = initVal[kk];
			}
			k = 1;
			while (k < maxIteration) {
				Arrays.fill(summ, 0.0);
				for (j = 0; j < i; j++) {
					rhsMat.mul(x, y);
					temp = MVO.dot_product(eigenVectors[j], y);
					for (int jj = 0; jj < n; jj++) {
						summ[jj] += eigenVectors[j][jj] * temp;
					}
				}
				for (int ii = 0; ii < n; ii++) {
					x[ii] -= summ[ii];
				}
				rhsMat.mul(x, y);
				double xMx = Math.sqrt(MVO.dot_product(x, y));
				for (int ii = 0; ii < n; ii++) {
					x[ii] /= xMx;
				}
				rhsMat.mul(x, y);
				//solver.PCGSSOR(y, xx, 1.5, 0);
				xx = solver.Solve_LU(y, 0);
				System.out.println("Find\t" + (i + 1) + "\t of\t" + numOfEV + "\t eigenvalues!");
				KK = MVO.dot_product(xx, y);
				rhsMat.mul(xx, y);
				MM = MVO.dot_product(xx, y);
				rr = KK / MM;
				if (k > 1 && Math.abs((rr - r) / rr) < this.eps) {
					break;
				}
				k++;
				r = rr;
				for (int kk = 0; kk < n; kk++) {
					x[kk] = xx[kk];
				}

			}
			System.out.println((i + 1) + "  of  " + numOfEV + "  finished!");
			// write(*, *) rr - alpha
			eigenValues[i] = rr - offsetVal;
			for (int kk = 0; kk < n; kk++) {
				eigenVectors[i][kk] = xx[kk] / Math.sqrt(MM);
			}
			System.gc();
		}
	}
}
